home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / dsterf.f < prev    next >
Text File  |  1997-06-25  |  10KB  |  382 lines
  1.       SUBROUTINE DSTERF( N, D, E, INFO )
  2. *
  3. *  -- LAPACK routine (version 2.0) --
  4. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  5. *     Courant Institute, Argonne National Lab, and Rice University
  6. *     September 30, 1994
  7. *
  8. *     .. Scalar Arguments ..
  9.       INTEGER            INFO, N
  10. *     ..
  11. *     .. Array Arguments ..
  12.       DOUBLE PRECISION   D( * ), E( * )
  13. *     ..
  14. *
  15. *  Purpose
  16. *  =======
  17. *
  18. *  DSTERF computes all eigenvalues of a symmetric tridiagonal matrix
  19. *  using the Pal-Walker-Kahan variant of the QL or QR algorithm.
  20. *
  21. *  Arguments
  22. *  =========
  23. *
  24. *  N       (input) INTEGER
  25. *          The order of the matrix.  N >= 0.
  26. *
  27. *  D       (input/output) DOUBLE PRECISION array, dimension (N)
  28. *          On entry, the n diagonal elements of the tridiagonal matrix.
  29. *          On exit, if INFO = 0, the eigenvalues in ascending order.
  30. *
  31. *  E       (input/output) DOUBLE PRECISION array, dimension (N-1)
  32. *          On entry, the (n-1) subdiagonal elements of the tridiagonal
  33. *          matrix.
  34. *          On exit, E has been destroyed.
  35. *
  36. *  INFO    (output) INTEGER
  37. *          = 0:  successful exit
  38. *          < 0:  if INFO = -i, the i-th argument had an illegal value
  39. *          > 0:  the algorithm failed to find all of the eigenvalues in
  40. *                a total of 30*N iterations; if INFO = i, then i
  41. *                elements of E have not converged to zero.
  42. *
  43. *  =====================================================================
  44. *
  45. *     .. Parameters ..
  46.       DOUBLE PRECISION   ZERO, ONE, TWO, THREE
  47.       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
  48.      $                   THREE = 3.0D0 )
  49.       INTEGER            MAXIT
  50.       PARAMETER          ( MAXIT = 30 )
  51. *     ..
  52. *     .. Local Scalars ..
  53.       INTEGER            I, ISCALE, JTOT, L, L1, LEND, LENDM1, LENDP1,
  54.      $                   LENDSV, LM1, LSV, M, MM1, NM1, NMAXIT
  55.       DOUBLE PRECISION   ALPHA, ANORM, BB, C, EPS, EPS2, GAMMA, OLDC,
  56.      $                   OLDGAM, P, R, RT1, RT2, RTE, S, SAFMAX, SAFMIN,
  57.      $                   SIGMA, SSFMAX, SSFMIN, TST
  58. *     ..
  59. *     .. External Functions ..
  60.       DOUBLE PRECISION   DLAMCH, DLANST, DLAPY2
  61.       EXTERNAL           DLAMCH, DLANST, DLAPY2
  62. *     ..
  63. *     .. External Subroutines ..
  64.       EXTERNAL           DLAE2, DLASCL, DLASRT, XERBLA
  65. *     ..
  66. *     .. Intrinsic Functions ..
  67.       INTRINSIC          ABS, SIGN, SQRT
  68. *     ..
  69. *     .. Executable Statements ..
  70. *
  71. *     Test the input parameters.
  72. *
  73.       INFO = 0
  74. *
  75. *     Quick return if possible
  76. *
  77.       IF( N.LT.0 ) THEN
  78.          INFO = -1
  79.          CALL XERBLA( 'DSTERF', -INFO )
  80.          RETURN
  81.       END IF
  82.       IF( N.LE.1 )
  83.      $   RETURN
  84. *
  85. *     Determine the unit roundoff for this environment.
  86. *
  87.       EPS = DLAMCH( 'E' )
  88.       EPS2 = EPS**2
  89.       SAFMIN = DLAMCH( 'S' )
  90.       SAFMAX = ONE / SAFMIN
  91.       SSFMAX = SQRT( SAFMAX ) / THREE
  92.       SSFMIN = SQRT( SAFMIN ) / EPS2
  93. *
  94. *     Compute the eigenvalues of the tridiagonal matrix.
  95. *
  96.       NMAXIT = N*MAXIT
  97.       SIGMA = ZERO
  98.       JTOT = 0
  99. *
  100. *     Determine where the matrix splits and choose QL or QR iteration
  101. *     for each block, according to whether top or bottom diagonal
  102. *     element is smaller.
  103. *
  104.       L1 = 1
  105.       NM1 = N - 1
  106. *
  107.    10 CONTINUE
  108.       IF( L1.GT.N )
  109.      $   GO TO 170
  110.       IF( L1.GT.1 )
  111.      $   E( L1-1 ) = ZERO
  112.       IF( L1.LE.NM1 ) THEN
  113.          DO 20 M = L1, NM1
  114.             TST = ABS( E( M ) )
  115.             IF( TST.EQ.ZERO )
  116.      $         GO TO 30
  117.             IF( TST.LE.( SQRT( ABS( D( M ) ) )*SQRT( ABS( D( M+
  118.      $          1 ) ) ) )*EPS ) THEN
  119.                E( M ) = ZERO
  120.                GO TO 30
  121.             END IF
  122.    20    CONTINUE
  123.       END IF
  124.       M = N
  125. *
  126.    30 CONTINUE
  127.       L = L1
  128.       LSV = L
  129.       LEND = M
  130.       LENDSV = LEND
  131.       L1 = M + 1
  132.       IF( LEND.EQ.L )
  133.      $   GO TO 10
  134. *
  135. *     Scale submatrix in rows and columns L to LEND
  136. *
  137.       ANORM = DLANST( 'I', LEND-L+1, D( L ), E( L ) )
  138.       ISCALE = 0
  139.       IF( ANORM.GT.SSFMAX ) THEN
  140.          ISCALE = 1
  141.          CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L+1, 1, D( L ), N,
  142.      $                INFO )
  143.          CALL DLASCL( 'G', 0, 0, ANORM, SSFMAX, LEND-L, 1, E( L ), N,
  144.      $                INFO )
  145.       ELSE IF( ANORM.LT.SSFMIN ) THEN
  146.          ISCALE = 2
  147.          CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L+1, 1, D( L ), N,
  148.      $                INFO )
  149.          CALL DLASCL( 'G', 0, 0, ANORM, SSFMIN, LEND-L, 1, E( L ), N,
  150.      $                INFO )
  151.       END IF
  152. *
  153.       DO 40 I = L, LEND - 1
  154.          E( I ) = E( I )**2
  155.    40 CONTINUE
  156. *
  157. *     Choose between QL and QR iteration
  158. *
  159.       IF( ABS( D( LEND ) ).LT.ABS( D( L ) ) ) THEN
  160.          LEND = LSV
  161.          L = LENDSV
  162.       END IF
  163. *
  164.       IF( LEND.GE.L ) THEN
  165. *
  166. *        QL Iteration
  167. *
  168. *        Look for small subdiagonal element.
  169. *
  170.    50    CONTINUE
  171.          IF( L.NE.LEND ) THEN
  172.             LENDM1 = LEND - 1
  173.             DO 60 M = L, LENDM1
  174.                TST = ABS( E( M ) )
  175.                IF( TST.LE.EPS2*ABS( D( M )*D( M+1 ) ) )
  176.      $            GO TO 70
  177.    60       CONTINUE
  178.          END IF
  179. *
  180.          M = LEND
  181. *
  182.    70    CONTINUE
  183.          IF( M.LT.LEND )
  184.      $      E( M ) = ZERO
  185.          P = D( L )
  186.          IF( M.EQ.L )
  187.      $      GO TO 90
  188. *
  189. *        If remaining matrix is 2 by 2, use DLAE2 to compute its
  190. *        eigenvalues.
  191. *
  192.          IF( M.EQ.L+1 ) THEN
  193.             RTE = SQRT( E( L ) )
  194.             CALL DLAE2( D( L ), RTE, D( L+1 ), RT1, RT2 )
  195.             D( L ) = RT1
  196.             D( L+1 ) = RT2
  197.             E( L ) = ZERO
  198.             L = L + 2
  199.             IF( L.LE.LEND )
  200.      $         GO TO 50
  201.             GO TO 150
  202.          END IF
  203. *
  204.          IF( JTOT.EQ.NMAXIT )
  205.      $      GO TO 150
  206.          JTOT = JTOT + 1
  207. *
  208. *        Form shift.
  209. *
  210.          RTE = SQRT( E( L ) )
  211.          SIGMA = ( D( L+1 )-P ) / ( TWO*RTE )
  212.          R = DLAPY2( SIGMA, ONE )
  213.          SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) )
  214. *
  215.          C = ONE
  216.          S = ZERO
  217.          GAMMA = D( M ) - SIGMA
  218.          P = GAMMA*GAMMA
  219. *
  220. *        Inner loop
  221. *
  222.          MM1 = M - 1
  223.          DO 80 I = MM1, L, -1
  224.             BB = E( I )
  225.             R = P + BB
  226.             IF( I.NE.M-1 )
  227.      $         E( I+1 ) = S*R
  228.             OLDC = C
  229.             C = P / R
  230.             S = BB / R
  231.             OLDGAM = GAMMA
  232.             ALPHA = D( I )
  233.             GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM
  234.             D( I+1 ) = OLDGAM + ( ALPHA-GAMMA )
  235.             IF( C.NE.ZERO ) THEN
  236.                P = ( GAMMA*GAMMA ) / C
  237.             ELSE
  238.                P = OLDC*BB
  239.             END IF
  240.    80    CONTINUE
  241. *
  242.          E( L ) = S*P
  243.          D( L ) = SIGMA + GAMMA
  244.          GO TO 50
  245. *
  246. *        Eigenvalue found.
  247. *
  248.    90    CONTINUE
  249.          D( L ) = P
  250. *
  251.          L = L + 1
  252.          IF( L.LE.LEND )
  253.      $      GO TO 50
  254.          GO TO 150
  255. *
  256.       ELSE
  257. *
  258. *        QR Iteration
  259. *
  260. *        Look for small superdiagonal element.
  261. *
  262.   100    CONTINUE
  263.          IF( L.NE.LEND ) THEN
  264.             LENDP1 = LEND + 1
  265.             DO 110 M = L, LENDP1, -1
  266.                TST = ABS( E( M-1 ) )
  267.                IF( TST.LE.EPS2*ABS( D( M )*D( M-1 ) ) )
  268.      $            GO TO 120
  269.   110       CONTINUE
  270.          END IF
  271. *
  272.          M = LEND
  273. *
  274.   120    CONTINUE
  275.          IF( M.GT.LEND )
  276.      $      E( M-1 ) = ZERO
  277.          P = D( L )
  278.          IF( M.EQ.L )
  279.      $      GO TO 140
  280. *
  281. *        If remaining matrix is 2 by 2, use DLAE2 to compute its
  282. *        eigenvalues.
  283. *
  284.          IF( M.EQ.L-1 ) THEN
  285.             RTE = SQRT( E( L-1 ) )
  286.             CALL DLAE2( D( L ), RTE, D( L-1 ), RT1, RT2 )
  287.             D( L ) = RT1
  288.             D( L-1 ) = RT2
  289.             E( L-1 ) = ZERO
  290.             L = L - 2
  291.             IF( L.GE.LEND )
  292.      $         GO TO 100
  293.             GO TO 150
  294.          END IF
  295. *
  296.          IF( JTOT.EQ.NMAXIT )
  297.      $      GO TO 150
  298.          JTOT = JTOT + 1
  299. *
  300. *        Form shift.
  301. *
  302.          RTE = SQRT( E( L-1 ) )
  303.          SIGMA = ( D( L-1 )-P ) / ( TWO*RTE )
  304.          R = DLAPY2( SIGMA, ONE )
  305.          SIGMA = P - ( RTE / ( SIGMA+SIGN( R, SIGMA ) ) )
  306. *
  307.          C = ONE
  308.          S = ZERO
  309.          GAMMA = D( M ) - SIGMA
  310.          P = GAMMA*GAMMA
  311. *
  312. *        Inner loop
  313. *
  314.          LM1 = L - 1
  315.          DO 130 I = M, LM1
  316.             BB = E( I )
  317.             R = P + BB
  318.             IF( I.NE.M )
  319.      $         E( I-1 ) = S*R
  320.             OLDC = C
  321.             C = P / R
  322.             S = BB / R
  323.             OLDGAM = GAMMA
  324.             ALPHA = D( I+1 )
  325.             GAMMA = C*( ALPHA-SIGMA ) - S*OLDGAM
  326.             D( I ) = OLDGAM + ( ALPHA-GAMMA )
  327.             IF( C.NE.ZERO ) THEN
  328.                P = ( GAMMA*GAMMA ) / C
  329.             ELSE
  330.                P = OLDC*BB
  331.             END IF
  332.   130    CONTINUE
  333. *
  334.          E( LM1 ) = S*P
  335.          D( L ) = SIGMA + GAMMA
  336.          GO TO 100
  337. *
  338. *        Eigenvalue found.
  339. *
  340.   140    CONTINUE
  341.          D( L ) = P
  342. *
  343.          L = L - 1
  344.          IF( L.GE.LEND )
  345.      $      GO TO 100
  346.          GO TO 150
  347. *
  348.       END IF
  349. *
  350. *     Undo scaling if necessary
  351. *
  352.   150 CONTINUE
  353.       IF( ISCALE.EQ.1 )
  354.      $   CALL DLASCL( 'G', 0, 0, SSFMAX, ANORM, LENDSV-LSV+1, 1,
  355.      $                D( LSV ), N, INFO )
  356.       IF( ISCALE.EQ.2 )
  357.      $   CALL DLASCL( 'G', 0, 0, SSFMIN, ANORM, LENDSV-LSV+1, 1,
  358.      $                D( LSV ), N, INFO )
  359. *
  360. *     Check for no convergence to an eigenvalue after a total
  361. *     of N*MAXIT iterations.
  362. *
  363.       IF( JTOT.EQ.NMAXIT ) THEN
  364.          DO 160 I = 1, N - 1
  365.             IF( E( I ).NE.ZERO )
  366.      $         INFO = INFO + 1
  367.   160    CONTINUE
  368.          RETURN
  369.       END IF
  370.       GO TO 10
  371. *
  372. *     Sort eigenvalues in increasing order.
  373. *
  374.   170 CONTINUE
  375.       CALL DLASRT( 'I', N, D, INFO )
  376. *
  377.       RETURN
  378. *
  379. *     End of DSTERF
  380. *
  381.       END
  382.